{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "The [HHL algorithm](https://en.wikipedia.org/wiki/Quantum_algorithm_for_linear_systems_of_equations) underlies many quantum machine learning protocols, but it is a highly nontrivial algorithm with lots of conditions. In this notebook, we implement the algorithm to gain a better understanding of how it works and when it works efficiently. The notebook is derived from the [computational appendix](https://gitlab.com/apozas/bayesian-dl-quantum) of the paper [Bayesian Deep Learning on a Quantum Computer](https://arxiv.org/abs/1806.11463). We restrict our attention to inverting a $2\\times 2$ matrix, following Pan *et al*.'s implementation [[1](#1)] of the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:39.391070Z",
     "start_time": "2018-11-19T20:10:38.023338Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/pwittek/.anaconda3/envs/qiskit/lib/python3.7/site-packages/marshmallow/schema.py:364: ChangedInMarshmallow3Warning: strict=False is not recommended. In marshmallow 3.0, schemas will always be strict. See https://marshmallow.readthedocs.io/en/latest/upgrading.html#schemas-are-always-strict\n",
      "  ChangedInMarshmallow3Warning\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister\n",
    "from qiskit import execute\n",
    "from qiskit import BasicAer\n",
    "π = np.pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setting up the problem\n",
    "\n",
    "We will solve the equation $Ax=b$ with $A = \\frac{1}{2}\\begin{bmatrix}3 & 1 \\\\1 & 3 \\\\ \\end{bmatrix}$ and $b =\\begin{bmatrix} 1 \\\\ 0 \\\\ \\end{bmatrix}$. We will encode the $A$ matrix as a Hamiltonian and $b$ in a register. With ancillas, we will need a total of five qubits and one classical register for post-selection. We add an extra qubit and extra classical register to create a swap test to compare our result to the ideal state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:39.451554Z",
     "start_time": "2018-11-19T20:10:39.393922Z"
    }
   },
   "outputs": [],
   "source": [
    "q = QuantumRegister(6)\n",
    "c = ClassicalRegister(2)\n",
    "hhl = QuantumCircuit(q, c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vector $b$ can be encoded as $\\left|b\\right\\rangle = \\sum_{i=0}^N b_i\\left|i\\right\\rangle = \\left|0\\right\\rangle$, so no explicit state preparation circuit is needed for this case (this will not be true in general).\n",
    "\n",
    "# Quantum phase estimation\n",
    "\n",
    "The next step is to encode the eigenvalues of the matrix $A$ in an additional register. This is done via quantum phase estimation of the evolution described by the Hamiltonian $A$ during some time $t_0$, $\\exp(\\imath A t_0)$. The protocol has three steps.\n",
    "\n",
    "First we prepare the ancilla state $\\left|\\psi_0\\right\\rangle=\\sum_{\\tau=0}^{T-1}\\left|\\tau\\right\\rangle$. Why this state? It will control the time evolution: it is like a clock, turning on evolution for a certain amount of time. The original HHL algorithm suggests a weighted superposition of all states $\\tau$ that minimizes errors in following steps in the algorithm. However, for our implementation, a uniform superposition already gives good results.\n",
    "\n",
    "Our goal is to create a superposition of $A$ as a Hamiltonian applied for different durations. Since the eigenvalues are always situated on the complex unit circle, these differently evolved components in the superposition help reveal the eigenstructure. So we apply the conditional Hamiltonian evolution $\\sum_{\\tau=0}^{T-1}\\left|\\tau\\right\\rangle\\left\\langle\\tau\\right|\\otimes e^{\\imath A\\tau t_0/T}$ on $\\left|\\psi_0\\right\\rangle\\otimes\\left|b\\right\\rangle$. This operation evolves the state $\\left|b\\right\\rangle$ according to the Hamiltonian $A$ for the time $\\tau$ determined by the state $\\left|\\psi_0\\right\\rangle$. Given that in $\\left|\\psi_0\\right\\rangle$ we have a superposition of all possible time steps between $0$ and $T$, we will end up with a superposition of all possible evolutions, and a suitable choice of number of timesteps $T$ and total evolution time $t_0$ allow to encode binary representations of the eigenvalues.\n",
    "\n",
    "As a final step, we apply an inverse Fourier transformation that writes the phases (that, recall, encode the eigenvalues of $A$) into new registers.\n",
    "\n",
    "The total circuit for phase estimation is the following:\n",
    "\n",
    "![Quantum phase estimation in the quantum matrix inversion algorithm](figures/qpe_for_hhl.svg)\n",
    "\n",
    "\n",
    "In our $2\\times 2$ case, the circuit is massively simplified. Given that the matrix $A$ has eigenvalues that are powers of $2$, we can choose $T=4$, $t_0=2\\pi$ to obtain exact results with just two controlled evolutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:39.504410Z",
     "start_time": "2018-11-19T20:10:39.457938Z"
    }
   },
   "outputs": [],
   "source": [
    "# Superposition\n",
    "hhl.h(q[1])\n",
    "hhl.h(q[2])\n",
    "# Controlled-U0\n",
    "hhl.cu3(-π / 2, -π / 2, π / 2, q[2], q[3])\n",
    "hhl.cu1(3 * π / 4, q[2], q[3])\n",
    "hhl.cx(q[2], q[3])\n",
    "hhl.cu1(3 * π / 4, q[2], q[3])\n",
    "hhl.cx(q[2], q[3])\n",
    "# Controlled-U1\n",
    "hhl.cx(q[1], q[3]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We apply quantum inverse Fourier transformation to write the phase to a register:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:39.540991Z",
     "start_time": "2018-11-19T20:10:39.522622Z"
    }
   },
   "outputs": [],
   "source": [
    "hhl.swap(q[1], q[2])\n",
    "hhl.h(q[2])\n",
    "hhl.cu1(-π / 2, q[1], q[2])\n",
    "hhl.h(q[1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The state of the system after this decomposition is approximately $\\sum _{j{\\mathop {=}}1}^{N}\\beta _{j}\\left|u_{j}\\right\\rangle \\left|\\lambda_{j}\\right\\rangle$, where $\\left|b\\right\\rangle=\\sum _{j{\\mathop {=}}1}^{N}\\beta _{j}\\left|u_{j}\\right\\rangle$ is the encoding of the vector $b$ in the eigenbasis of $A$. Now, there is an often overlooked step that performs bit operations on $\\left|\\lambda_{j}\\right\\rangle$ to actually invert it.\n",
    "\n",
    "In our case, the inversion of the eigenvalues is easy. The eigenvalues of $A$ are $\\lambda_1=2=10_2$ and $\\lambda_2=1=01_2$, and their reciprocals are $\\lambda_1^{-1}=1/2$ and $\\lambda_2^{-1}=1$. Noting that $2\\lambda_1^{-1}=01_2$ and $2\\lambda_2^{-1}=10_2$, a swap gate is enough to obtain the state $\\sum _{j{\\mathop {=}}1}^{N}\\beta _{j}\\left|u_{j}\\right\\rangle \\left|2\\lambda _{j}^{-1}\\right\\rangle$, that encodes the reciprocals of the eigenvalues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:39.557636Z",
     "start_time": "2018-11-19T20:10:39.543737Z"
    }
   },
   "outputs": [],
   "source": [
    "hhl.swap(q[1], q[2]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conditional rotation of ancilla\n",
    "\n",
    "Next, we perform a conditional rotation to encode the information of the reciprocals of the eigenvalues in the amplitudes of a state, on which we will later post-select. The state we would like to get is $\\sum _{j{\\mathop {=}}1}^{N}\\beta _{j}\\left|u_{j}\\right\\rangle\\left|2\\lambda _{j}^{-1}\\right\\rangle \\left(\\sqrt{1-\\frac{C^2}{\\lambda_j^2}}\\left|0\\right\\rangle+\\frac{C}{\\lambda_j}\\left|1\\right\\rangle \\right)$. This is achieved by controlled rotations in the same spirit of the conditional Hamiltonian evolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:39.579455Z",
     "start_time": "2018-11-19T20:10:39.560300Z"
    }
   },
   "outputs": [],
   "source": [
    "hhl.cu3(0.392699, 0, 0, q[1], q[0])  # Controlled-RY0\n",
    "hhl.cu3(0.19634955, 0, 0, q[2], q[0]);  # Controlled-RY1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Uncomputing the eigenvalue register\n",
    "\n",
    "A necessary step when performing quantum computations is to uncompute all operations except those that store the information that we want to obtain from the algorithm in the final registers. We need to do this in case the registers are entangled, which would affect the results.\n",
    "\n",
    "In our case, we must uncompute the phase estimation protocol. After the uncomputation, the state should be $\\sum_{j=1}^N\\beta_j\\left|u_j\\right\\rangle\\left|0\\right\\rangle\\left(\\sqrt{1-\\frac{C^2}{\\lambda_j^2}}\\left|0\\right\\rangle+\\frac{C}{\\lambda_j}\\left|1\\right\\rangle \\right)$, so we can safely forget about the eigenvalue register."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:39.620961Z",
     "start_time": "2018-11-19T20:10:39.581605Z"
    }
   },
   "outputs": [],
   "source": [
    "hhl.swap(q[1], q[2])\n",
    "hhl.h(q[1])\n",
    "hhl.cu1(π / 2, q[1], q[2])  # Inverse(Dagger(Controlled-S))\n",
    "hhl.h(q[2])\n",
    "hhl.swap(q[2], q[1])\n",
    "# Inverse(Controlled-U1)\n",
    "hhl.cx(q[1], q[3])\n",
    "# Inverse(Controlled-U0)\n",
    "hhl.cx(q[2], q[3])\n",
    "hhl.cu1(-3 * π / 4, q[2], q[3])\n",
    "hhl.cx(q[2], q[3])\n",
    "hhl.cu1(-3 * π / 4, q[2], q[3])\n",
    "hhl.cu3(-π / 2, π / 2, -π / 2, q[2], q[3])\n",
    "# End of Inverse(Controlled-U0)\n",
    "hhl.h(q[2])\n",
    "hhl.h(q[1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rejection sampling on the ancilla register and a swap test\n",
    "\n",
    "The state $\\left|x\\right\\rangle=A^{-1}\\left|b\\right\\rangle\\propto\\sum_j \\beta_j\\lambda_j^{-1}\\left|u_j\\right\\rangle$ that contains information about the solution to $Ax=b$ is that obtained when measuring $1$ on the ancilla state. We perform the post-selection by projecting onto the desired $\\left|1\\right\\rangle$. To check that the solution is the expected one, we prepare the correct output state manually to perform a swap test with the outcome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:39.654681Z",
     "start_time": "2018-11-19T20:10:39.623580Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Target state preparation\n",
    "hhl.rz(-π, q[4])\n",
    "hhl.u1(π, q[4])\n",
    "hhl.h(q[4])\n",
    "hhl.ry(-0.9311623288419387, q[4])\n",
    "hhl.rz(π, q[4])\n",
    "# Swap test\n",
    "hhl.h(q[5])\n",
    "hhl.cx(q[4], q[3])\n",
    "hhl.ccx(q[5], q[3], q[4])\n",
    "hhl.cx(q[4], q[3])\n",
    "hhl.h(q[5])\n",
    "\n",
    "hhl.barrier(q)\n",
    "hhl.measure(q[0], c[0])\n",
    "hhl.measure(q[5], c[1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note: it is a good exercise to check that the right result is given by the state $\\left|x\\right\\rangle=0.949\\left|0\\right\\rangle + 0.314\\left|1\\right\\rangle$, which is the state we prepare above.*\n",
    "\n",
    "There are two measurements performed, one of the ancilla register (for doing the post-selection) and another one that gives the result of the swap test. To calculate success probabilities, let us define some helper functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:39.671821Z",
     "start_time": "2018-11-19T20:10:39.656877Z"
    }
   },
   "outputs": [],
   "source": [
    "def get_psuccess(counts):\n",
    "    '''Compute the success probability of the HHL protocol from the statistics\n",
    "\n",
    "    :return: (float) The success probability.\n",
    "    '''\n",
    "    try:\n",
    "        succ_rotation_fail_swap = counts['11']\n",
    "    except KeyError:\n",
    "        succ_rotation_fail_swap = 0\n",
    "    try:\n",
    "        succ_rotation_succ_swap = counts['01']\n",
    "    except KeyError:\n",
    "        succ_rotation_succ_swap = 0\n",
    "    succ_rotation = succ_rotation_succ_swap + succ_rotation_fail_swap\n",
    "    try:\n",
    "        prob_swap_test_success = succ_rotation_succ_swap / succ_rotation\n",
    "    except ZeroDivisionError:\n",
    "        prob_swap_test_success = 0\n",
    "    return prob_swap_test_success"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we run the circuit on the simulator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:42.138218Z",
     "start_time": "2018-11-19T20:10:39.673729Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n"
     ]
    }
   ],
   "source": [
    "backend = BasicAer.get_backend('qasm_simulator')\n",
    "job = execute(hhl, backend, shots=100)\n",
    "result = job.result()\n",
    "counts = result.get_counts(hhl)\n",
    "print(get_psuccess(counts))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running on the actual QPU would yield a much poorer result due to imprecisions in the applications of the gates and noise caused by the environment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "[1] J. Pan, Y. Cao, X. Yao, Z. Li, C. Ju, H. Chen, X. Peng, S. Kais, and J. Du. (2014). [Experimental realization of quantum algorithm for solving linear systems of equations](https://arxiv.org/abs/1302.1946). *Physical Review Letters* 89:022313. <a id='1'></a>\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
